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In this work, optimal shaft shapes for flow 
in the annular space between a rotating 
shaft with axially-periodic radius and a 
fixed coaxial outer circular cylinder, are 
investigated. Axisymmetric steady flows 
in this geometry are determined by solving 
the full Navier-Stokes equations in the 
actual domain. A measure of the flow 
field, a weighted convex combination 
of the volume averaged square of the 
/,2-iiorm of the velocity and vorticity 
vectors, is employed. It has been 
demonstrated that boundary shape can 
be used to influence the characteristics 
of the flow field, such as its velocity 
component distribution, kinetic energy, 
or even vorticity. This ability to influence 
flow fields through boundary shape may 
be employed to improve microfluidic 
mixing or, possibly, to minimize shear in 
biological applications. 
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1. Introduction 

Flow control and optimization has been an experimen- 
tal science for well over 100 years, with an early exam- 
ple being the work by Prandtl [1] in the early 1900s. The 
bulk of the early work on flow control and optimization 
fall into two categories: attempts to control fluid motion 
without the use of sophisticated fluid models involving 
partial differential equations (e.g., design of heating and 
cooling systems in a building) and attempts to control 
fluid motion without the use of sophisticated optimiza- 
tion techniques (e.g., boundary layer control), but using 
sophisticated fluid models. More recently, because of a 
better understanding of the equations governing fluid 
flow, more robust and better optimization algorithms 
have been employed, together with sophisticated fluid 
models to solve very complicated flow control and flow 
design problems, (see for example [2-3] ). 



Flow in or through a circular annulus driven by rota- 
tion of the inner cylinder has been of interest since the 
work of Taylor [4]. One way to influence this flow is 
to change the geometry of either the inner or outer cylin- 
der. Rotation of a shaft with axially- varying radius gives 
rise to a radial velocity component at any non-zero rota- 
tion rate, which in turn provides convective mixing in the 
radial direction. It has been shown that for the case of 
axisymmetric, axially-periodic variation of the inner 
cylinder radius, one can achieve significant modification 
of the flow [5]. By contrast, at rotation rates below the 
onset of Taylor vortices, a constant-radius shaft provides 
no radial mixing beyond diffusion. This, and the absence 
of sharp edges, may make this approach attractive for 
applications where cell damage must be avoided; thus 
controlling this flow is of interest. 

Selecting an objective function to optimize is an 
application and problem specific challenge. The 
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governing Navier-Stokes equations are highly nonlin- 
ear partial differential equations (quadratic in the veloc- 
ity), making simulation a computationally demanding 
task. Moreover, provisions must be made for changes in 
the shape of the boundary. This can be dealt with by 
using a grid that includes the entire domain (i.e., grid 
the solid and fluid portions of the domain) and a grid 
spacing that does not change from one geometry to 
another [2]. This method could be inefficient in that 
grid points in the solid portion of the domain would 
either need to be removed by algebraic manipulation or 
have their velocity components set to zero. A fictitious 
domain method could be employed [2]. 

In simulation-based optimization settings, a comput- 
er simulation must be run repeatedly, in order to com- 
pute the various quantities needed by the optimization 
algorithm. The resulting simulation output must then be 
postprocessed in order to evaluate the cost fimction. 
The cost of computing the cost fimction is small when 
compared with that of obtaining derivatives for gradi- 
ent-based methods, but in many cases these derivatives 
are not accurate, not available, or simply too costly to 
calculate [6]. Although automatic differentiation tools 
have been shown to be effective for some types of 
simulation-based problems, they are not universally 
applicable [7]. Results of complex calculations may fail 
to have the level of precision necessary for a reliable 
finite-difference approximation to the gradient, elimi- 
nating Newton or quasi-Newton methods as possible 
optimization algorithms. In general, features such as 
adaptive algorithms, stopping tests in iterative schemes 
inside the simulation, and the inevitable effect of float- 
ing point arithmetic are among the culprits that cause 
smooth problems to appear nonsmooth [6]. For these 
reasons, and the fact that one can parameterize a large 
class of shapes with a small number of parameters (in 
this case no more than five), a direct search method can 
be employed to efficiently and effectively find optimal 
shaft shapes. 

The paper is organized as follows: the fluid dynam- 
ics problem is reviewed in §2, followed by a descrip- 
tion of the optimization approach in §3. Results and 
discussion are presented in §4, followed by conclusions 
in §5. 

2. Flow 

This section provides a brief review of the formula- 
tion and the numerical method used in computing 
axisymmetric solutions (i.e., solutions that are invariant 



along the azimuthal coordinate) of the incompressible 
Navier-Stokes equations for flow between a fixed outer 
circular cylinder and a rotating coaxial inner shaft with 
either rectangular or sinusoidal shape. 

Work on extended (i.e., periodic or random) radius 
variation includes experimental work [8-19], while 
related work on the effect of localized radius variation 
has been reported in [20-22]. The most detailed of the 
experimental investigations are those of [8] and [9] (see 
[5] for a detailed review of previous work). On the 
computational side, the most detailed calculations are 
those of [5], in which bifiircation behavior results as a 
fimction of flow rate were presented for a radius ratio 
of 7] = R^/Ro = 0.5 (where R^ and Rq are defined 
below) and several values of the radial amplitude. 

2.1 Formulation 

In this section, we define the governing equations 
without specifying the exact shape of the boundary. In 
general, the geometry of interest is that of a fixed outer 
circular cylinder and an inner coaxial rotating shaft 
with topography. Axisymmetric steady flows in this 
geometry are determined by solving the fiill Navier- 
Stokes equations in the actual domain. The governing 
equations and boundary conditions are nondimension- 
alized by scaling the dimensional coordinates R and Z 
with the outer radius, Rg, and the axial wavelength, L, 
respectively while the velocity and pressure are nondi- 
mensionalized using Q^R^ and pO^Rl respectively; 
where p is the fluid density, Q^ is the angular velocity 
of the inner shaft, and Ri is an appropriate measure of 
the inner radius (e.g., the mean inner shaft radius). 
Radius R^ will be defined more precisely, later in this 
section when the shaft cross-sectional shapes are 
introduced. 

In an inertial reference frame, the dimensionless 
continuity and momentum equations are 



dv V ^dv 

dr r dz 



(la) 



^^-^+i(v^H-4|(lb) 



dr dz r dr Ta 



dv^ dv^ v^v^ 1 ( v^) . s 



dv^ o dv^ ndp 1 ^2 



(Id) 
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where 



r dr dr dz 



(2) 



I3 = Rq/L is the length scale ratio, and Ta = (0^ R^ Rq)/v 
is the Taylor number, where v is the kinematic viscosi- 
ty of the fluid. 

We enforce no-slip conditions on the rigid imperme- 
able inner shaft and outer circular cylinder 

v^ir.{z),z)=0 and v^r^^z) =0 (3a,b) 



v^ir.(z%z) = and l>^(ro,z)=0 



(3c,d) 



^^(n(^\^) = f'i(^)^l and v^(r^,z)=0, (3e,f) 
and the flow is assumed to be axially periodic 

v(r,0) = v(r,l) (4a) 

3v(r,0) _ 3v(r,l) (4^^>l 



dz dz 

pir,0) = pir,l). 



(4c) 



hi this work two shaft shapes are considered; an 
azimuthal groove of rectangular cross section (Fig. la) 
and an azimuthal groove of sinusoidal cross section 
(Fig. lb). In a plane passing through the symmetry axis of 
the outer cylinder, Figs, la and lb show the two shaft geo- 
metries considered. In Fig. la, (^=(7?o-^i)/(^o-^i) is 
the clearance ratio, and T= IIL is the groove width ratio, 
where R^ is the radius of the flight and / is the groove 
width. This geometry is fully defined by setting ^8, r], f, 
and T. Figure lb shows the sinusoidal inner shaft geome- 
try with the inner radius varying according to 



r.{z) = ri + e^m{2nz). 



(5) 



where e is the amplitude of the radial modulation. This 
geometry is fully defined by setting ^8, r], and £. 

Note that the precise value of Ri depends on the 
geometry considered. For the rectangular groove case 
the radial length scale was chosen to be the hub radius 
(i.e., Ri = Rd, while for the sinusoidal groove case the 
mean inner radius (i.e., Ri = Rd is used. Thus, the pre- 
cise definitions of rj and Ta also depend on which 
geometry one is considering. 
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Fig. 1. Shaft geometry and computational domain, (a) azimuthal 
groove of rectangular shape, (b) azimuthal groove of sinusoidal 
shape. 



2.2 Numerical Method 

The momentum and continuity equations are solved 
numerically using a finite-element method, employing 
isoparametric quadrilateral elements with quadratic 
velocity and discontinuous linear pressure interpolation 
(cf [23-25]); a consistent penalty method to satisfy 
incompressibility Eq. (la) and a 5x5 Gaussian quadra- 
ture rule for numerical integration. A structured ortho- 
gonal mesh with the capability to adjust mesh density 
throughout the computational domain provides flexibil- 
ity to locally add more points to resolve relevant flow 
structures. 
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Discretization of Eq. (1) leads to a quadratically non- 
linear algebraic system of the form 



f(v) = 0, 



(6) 



We define the following objective function 



(8) 



which is solved by Newton iteration. If needed, con- 
verged solutions at one spatial resolution (or lower Ta) 
can be interpolated to get an initial iterate for the next- 
higher resolution. The convergence criteria were that 
both ||R|| and ||v"^' - v"|| / ||v"|| be less than 10"^, where 
III! denotes an Z2-norm [25], R is the residual vector 
used by the Newton iteration, and v is the velocity 
vector. A grid convergence study has been performed 
for various geometric parameter combinations span- 
ning the ranges considered in the work. Grid spacing 
gives well converged solutions (i.e., changing less than 
1 % with increased spatial resolution) for any para- 
meter combination considered. Results given in this 
paper are converged to the number of significant 
figures shown. See [5] for convergence results. 

3. Optimization 

Scalar measures of the flow and an objective func- 
tion are defined in this section. The problem of finding 
optimal shapes is considered in two stages; simulation 
of the fluid dynamics and parameter optimization. A 
finite element method, as described in the previous 
section, can be used to perform the simulation computa- 
tions, while a direct search software package, APPSPACK 
[26], is employed to estimate optimal flow parameters. 

3.1 Formulation 

Consider the volume averaged Z2-norm squared of 
the velocity (y^) and vorticity vectors (yj 



n 



r II ii2 

J II v|| r drdzdO 

V 

j r drdzdO 



(7a) 



r„=- 



w r drdzdO 



\r drdzdO 



(7b) 



where w is the vorticity vector. We note that these 
measures, for example, could be used to investigate 
mixing using only the velocity field, thus they have not 
been arbitrarily chosen. 



where J= ayy + {\ - a)/^^ is a convex combination of 
Eqs. (7a) and (7b), a is a constant taking on values 
between zero and one, J^ is a desired value of J, and F 
is a measure of the topography. We note that F differs 
between shaft shapes, with it being the dimensionless 
flight height, (1 - vi){\ - Q, for the rectangular groove 
and the dimensionless amplitude, £, for the sinusoidal 
groove. Minimumization of Eq. (8) over the geometric 
parameter space allows one to find the shaft shape that 
gives a value of J closest to the desired value, while 
also trying to keep the topography to a minimum. We 
note, for the sinusoidal groove case, that there exists a 
relationship between 77 and £ that must be taken into 
consideration. This relationship takes the form of a con- 
straint on £ 



0<£ <min(7?,l-7?). 

0<T7<1 



(9) 



Because the software package APPSPACK is 
employed, which handles only bound constraints, the 
constraint Eq. (9) must be handled by other means. In 
this work, constraints are handled through a special 
penalty function in the same way as employed by [2]. 
The penalty function penalizes infeasibility but returns 
a smaller function value for points that simultaneously 
reduce the objective function and infeasibility. In this 
way, the algorithm steers iterates away from points that 
are not physically possible (e.g., the gap width goes 
negative). 

4. Results and Discussion 

In this work, the simplest axisymmetric flow possi- 
ble is investigated and optimal shaft shapes are calcu- 
lated. For a smooth shaft, the location of the first bi- 
furcation from purely azimuthal flow to something 
more complex (e.g., for 77 = 0.5 the first bifurcation is 
to a steady and axisymmetric secondary flow and takes 
place at about 7a = 68.19) is strongly dependent on the 
value of the radius ratio. Significant modifications in 
solution topology have been shown to take place after 
the first bifurcation [5]. It is also known that the loca- 
tion of the first bifurcation decreases with increasing 
radius ratio [27-29]. Previous results [5] suggest that 
the location of the first bifurcation is significantly 
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affected by topography where for large amplitude 
cases results show that the location of the first bifurca- 
tion relocates to higher values of Ta (i.e., the modified 
Couette flow is much more stable than the smooth 
walled case). Computational results for 7a =10 are 
presented. This value was selected, in part, because it 
lies below the first bifurcation for any smooth shaft 
case for the combinations of geometric parameters 
considered here. This insures that solutions remain on a 
single solution branch, allowing for the optimal shaft 
shape to be solution branch independent (see [5] for 
solution branch details). 

In order to gage the effect (5 has on the flow field, we 
have defined an aspect ratio based on the gap and axial 
wavelength A = H/L, where H = {Rq - R^ for the rec- 
tangular azimuthal groove case and H=(Rq-R^ for 
the sinusoidal azimuthal groove case. This ratio can 
also be written in dimensionless form as A = j3 (1 - 77), 
where 7] is defined differently for each shaft shape, as 
discussed in §3. We note that figures are shown with 
respect to this aspect ratio. 

4.1 Rectangular Groove 

In this section, optimal shaft shapes for the rectangu- 
lar groove case are calculated. Geometric control 
parameters needed to fully define this geometry take on 
values in the following ranges 

j3e [0.5,2.5] and rye [0.2,0.8] (lOa, b) 
Ce [0.2,0.8] and re [0.2,0.8]. (lOc, d) 

For each a, F is minimized for both attainable and 
unattainable choices of J^, demonstrating that one can 
compute an optimal topography that generates this 
value without changing the rotation speed of the shaft. 
The results for a = 0.5 are qualitatively similar to those 
for a = 0; thus, results for a = 0.5 are not shown. 

a=0. For a=0, Figs. 2(a-d) show streamfunction 
and v^ (contour values range from to 1 . 1) contours for 
the set of geometric parameters that minimize and 
maximize J, and minimize F, given an attainable 
and unattainable J^, respectively. As seen in Fig. 2a 
(A = 0.35) for 13=0.5, rj = 0.294, C= 0.8, and T= 0.8, 
minimizing J (1 .62) leads to a shaft shape with minimal 
topography and small radius ratio. It is not surprising 
that the optimal shape would be one of large flight 
clearance and groove width, given that a measure of the 
vorticity field is being optimized. In the limit of zero 
flight height, the flow is unidirectional with only one 




Fig. 2. Streamfunction and Va contours for optimal shaft shapes 
for Ta= 10 and a= 0. a) Parameter set that minimizes /: j3= 0.5, 
77 = 0.29, C=0.8, and t=0.8. b) maximizes /: P = 2.5, 77 = 0.8, 
C=02, and t= 0.2. c) minimizes 7^ given an attainable J^: j3=0.59, 
77 = 0.2, ^= 0.2, and t= 0.5. d) minimizes F given an unattainable 
J^: P = 2.5, 77 = 0.8, C= 0.2, and t= 0.2. 



nonzero vorticity component at this Ta. The magnitude 
of azimuthal shear decreases with decreasing rj. As seen 
in Fig. 2b (A = 0.5) for j3 = 2.5, r] = 0.8, ^=0.2, and 
T= 0.2, maximizing J (473) leads again to a shaft shape 
with minimal topography. In contrast to the results 
shown in Fig. 2a when minimizing J, the topography 
maximizing J has a large radius ratio and small groove 
width. Because the azimuthal velocity component of the 
shaft goes as 13, R, having a large radius ratio and small- 
groove width maximizes the magnitude of this compo- 
nent over most of the axial length. 
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Topography can also fix the measure of the flow 
field to a desired value. This is of interest, for example, 
in mixing of biological agents, where excessive shear is 
known to cause damage. If one knows at what shear 
rate damage becomes significant, topography control 
can be used to avoid the maximum shear seen in the 
field, without changing the rotation rate. As seen in 
Fig. 2c (A = 0.48) for j3= 0.59, r] = 0.2, f = 0.2, t= 0.5, 
and an attainable J^ , minimizing F leads to a shaft shape 
that gives a value of /closest to the desired value while 
limiting changes in the topography to a minimum. For 
this case, F = 0.52 and F= 0.64. Thus, there is less than 
a 1 % difference between the J (200.45) that minimizes 



F and the desired value (J^ = 200). On the other hand, 
as shown in Fig. 2d {A = 0.5), minimizing F, given an 
unattainable value of J^ larger than possible at these 
rotation rates, leads to a shaft shape ((5 = 2.5, 7] = 0.8, 
^=0.2, T=0.2) that is identical to that found when 
maximizing J. 

a=l. For a= 1, Figs. 3(a-d) show streamfunction 
and t>0 (contour values range from to 1.1, with there 
being ten divisions) contours for the set of geometric 
parameters that minimize and maximize J, and mini- 
mize F, given an attainable and unattainable J^, respec- 
tively As seen in Fig. 3a (A = 0.4) for P = 0.5, rj = 0.2, 
^= 0.8, and T= 0.8, minimizing J once more leads to a 



a) 




b) 






L^) 



d) 




Fig. 3. Streamfunction and tJ^ contours for optimal shaft shapes for Ta = 10 and a= 1. a) Parameter set that minimizes J: p = 0.5, 
77 = 0.2, ^= 0.8, and t= 0.8. b) maximizes J: p= 2.5, 77= 0.2, C, = 0.2, and t= 0.33. c) minimizes F given an attainable J^\ p= 1.5, 
77 = 0.2, ^=0.43, and T= 0.35. d) minimizes F given an unattainable J^ : ^ =2.5, 77 = 0.2, f= 0.2, and t= 0.33. 
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shaft shape with minimal topography. The optimal shaft 
shape turns out to be one of large flight clearance and 
groove width, and small radius ratio. This is consistent 
with what one would expect when minimizing a meas- 
ure of the kinetic energy. As seen in Fig. 3b (A = 2) for 
j8=2.5, r] = 0.2, ^=0.2, and t=0.33, maximizing J 
leads to a shaft shape having significant topography 
(i.e., a deep thin groove). This differs from the a=0 
results, where the radius ratio was large. 

As shown in Fig. 3c (A= 1.2) for ^= 1.5, r] = 0.2, 
^= 0.43, T= 0.35, and an attainable J^, minimizing F 
leads to a shaft shape having medium flight clearance 
and small groove width. For this case, F=0.23 and 
r= 0.46. Also, for this case we have used the topo- 
graphy to exactly match the desired value of J^ = 3. We 
note that this value of ^ is significantly different than 
the values found for a = 0. The results (see Fig. 3d 
where A = 2) for an unattainable value of J^ are similar 
to those for a=0. 

4.2. Sinusoidal Groove 

In this section, optimal shaft shapes for the sinu- 
soidal groove case are calculated. Geometric para- 
meters needed to fully define this geometry are con- 
strained to take on values in the following ranges 

^e [0.5,2.5] and 7]^ [0.2,0.8] (lla,b) 

£e [0,0.5], (lie) 

where the upper bound on £ is taken to be 1/2 of the 
maximum value of the right-hand-side of Eq. (9). 



A minimum gap between the inner and outer cylinders 
of 0.05 is enforced. This is required for numerical rea- 
sons, because numerically it is desirable that a suffi- 
cient number of elements are present inside the gap to 
yield accurate solutions. The results for a =0.5 are 
again qualitatively similar to those for a = 0, as was the 
case for the rectangular groove case, and again are not 
shown. 

a= 0. For a=0, Figs. 4(a-b) show streamfunction 
and v^ (contour values range from to 1.1, with there 
being ten divisions) contours for the set of geometric 
parameters that maximize J and minimize F given an 
attainable J^, respectively. The optimal shaft shape that 
minimizes J has no topography. In the limit £ — > 0, the 
flow is unidirectional at this Ta and has only one nonze- 
ro vorticity component. For any nonzero value of £, 
there are three nonzero velocity components, and thus 
all three vorticity components are nonzero. As seen in 
Fig. 4a for P = 2.5, 7] = 0.70, and £ = 0.25, maximizing 
J leads to a shaft shape with significant topo- 
graphy. With A = 0.76, the mean gap is smaller than 
the axial wavelength. We note that as j8 increases 
while maintaining a constant mean gap width, the 
gradients in the flow also increase. This, along with 
the fact that w^ scales with P, causes y^ to increase 
as well. This optimal shaft shape is very similar to 
those found for the rectangular groove case, with 
j8 being the same, r] about 12 % less, and the minimum 
gap width chosen in each case. In contrast, the 
overall aspect ratio (A) is very different for these two 
cases. 



a) 




b) 




Fig. 4. Streamfunction and D^ contours for optimal sliaft sliapes for Ta= \0 and a = 0. a) Parameter set that maximizes J: p = 2.5, ?] = 0.70, and 
£= 0.25. b) minimizes F given an attainable J^'- P= 1.25, rf = 0.65, and £=0.25. 
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As seen in Fig. 4b for [5= 1.25, 7] = 0.65, 8=0.25, 
and an attainable J^ , minimizing F leads to a shaft shape 
that gives a value of /closest to the desired value, while 
keeping the topography to a minimum. For this case, 
F=0.17 and r= 0.25. As was the case for all the rec- 
tangular groove cases, there is less than a 1 % differ- 
ence between the J (25.09) that minimizes F and the 
desired value (J^ = 25). On the other hand, as was also 
the case for the rectangular groove case, minimizing F, 
given an unattainable value of J^ leads to a shaft shape 
(p = 2.5, 7] = 0.70, and £=0.25), identical to that found 
when maximizing J. 

a=l. For a= 1, Figs. 5(a-b) show streamfiinction 
and v^ (contour values range from to 1.1, with there 
being ten divisions) contours for the set of geometric 
parameters that maximize J and minimize F, given an 
attainable J^, respectively. As was the case for a = 0, the 
shaft shape that minimizes J has no topography. As seen 
in Fig. 5a for P = 2.5, rj = 0.35, and £ = 0.30, maximiz- 
ing J leads to a shaft shape with moderate topography 
compared to the a = case. As seen in Fig. 5a, A = 1 .63 
for this case, thus the mean gap width is much larger 
than the wavelength. 

As seen in Fig. 5b for p = 2.5, rj = O.S,e=0. 125, and 
an attainable J^, minimizing F leads to a shaft shape 
that gives a value of /closest to the desired value, while 
again keeping the topography to a minimum. For this 
case, F=O.OS and r= 0.125. There is about a 39% 
difference between the J (0.98) that minimizes F and 



a desired value (J^=0.6). This value is significantly 
larger than the value found for a = and any of the rec- 
tangular groove cases. For a= 1 , the range of J is much 
smaller than that found for a = 0, making it more diffi- 
cult to find the topography that gives the desired value. 
On the other hand, as was also the case for a = and all 
the rectangular groove cases, minimizing F, given an 
unattainable value of J^ leads to a shaft shape (^8= 2.5, 
?7 = 0.35, and £=0.30) identical to that found when 
maximizing J. 

5. Conclusions 

It has been demonstrated that boundary shape 
can be used to influence the characteristics of the 
flow field, such as its velocity component distribution, 
kinetic energy, or even vorticity The geometric 
parameters investigated in this work are those 
associated with mixing machinery and instruments 
being developed and used today. The work here 
provides a first step in designing an optimization tool 
for designing experiment and operation specific fluid 
mixing instruments. This ability to influence flow 
fields through boundary shape may be employed 
to improve microfludic mixing, maximize through 
flow while trying to improve low speed mixing, or 
possibly to minimize shear in biological applications 
(e.g., separations or downstream of artificial heart 
valves). 




b) 




Fig. 5. Streamfunction and Da contours for optimal shaft shapes for Ta= \0 and a = 1 . a) Parameter set that maximizes J: p = 2.5, rj = 0.35, and 
€= 0.30. b) mmimizes F given an attainable Jd- P = 2.5, 77 = 0.8, and € = 0. 125. 
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